The gut microbiome variability of a butterflyfish increases on severely degraded Caribbean reefs

Environmental degradation has the potential to alter key mutualisms that underlie the structure and function of ecological communities. How microbial communities associated with fishes vary across populations and in relation to habitat characteristics remains largely unknown despite their fundamental roles in host nutrition and immunity. We find significant differences in the gut microbiome composition of a facultative coral-feeding butterflyfish (Chaetodon capistratus) across Caribbean reefs that differ markedly in live coral cover (∼0–30%). Fish gut microbiomes were significantly more variable at degraded reefs, a pattern driven by changes in the relative abundance of the most common taxa potentially associated with stress. We also demonstrate that fish gut microbiomes on severely degraded reefs have a lower abundance of Endozoicomonas and a higher diversity of anaerobic fermentative bacteria, which may suggest a less coral dominated diet. The observed shifts in fish gut bacterial communities across the habitat gradient extend to a small set of potentially beneficial host associated bacteria (i.e., the core microbiome) suggesting essential fish-microbiome interactions may be vulnerable to severe coral degradation.

E nvironmental degradation associated with the Anthropocene is threatening the persistence of mutualistic relationships that are key to the stability of ecological functioning 1 . The increasingly severe degradation of coral reefs from both local and climatic stressors has led to novel habitat states with conspicuously altered fish and invertebrate communities [2][3][4][5] , making them a model system for studying ecological responses to environmental change. A potentially pervasive but largely overlooked response to habitat degradation is the change to host-associated microbiomes-the communities of bacteria, archaea, fungi, unicellular eukaryotes, protozoa, and viruses that live on internal and external surfaces of reef organisms. Host microbiomes potentially respond faster than their hosts to changing environmental conditions and can promote acclimatization processes as well as genetic adaptation [6][7][8][9][10] . Thus, microbial communities could play a key role in mediating a host's resilience and ability to adapt to environmental change. However, it remains to be explored whether mutualisms between fish hosts and gut microbiomes can shift to alternative beneficial relationships to provide a mechanism of resilience to habitat change, or whether the mutualism breaks down and simply reflects a cascading effect of degradation at all levels of ecological organization.
The importance of gut microbial communities in maintaining host health is well recognized in mammals and other vertebrates 11,12 , including a wealth of research into the importance of microbes in fish in aquaculture settings [13][14][15][16] . Fish harbor microbiomes that are unique from the microbial communities in their surrounding environment 17,18 . As the gut microbiome diversifies throughout the development of the fish host, a relatively stable gut microbiome is typically established within the first months of the fish's life 13 . These resident (autochthonous) microbes, which are consistently found associated with the fish population across space and time and potentially provide critical functions for the host, are referred to as the "core microbiome" 15,19,20 . In contrast, the numerous microbes occurring in the gastrointestinal tract after being ingested are transient (allochthonous) and may vary intraspecifically with the developmental stage and potentially include opportunistic pathogens. Because of their importance in maintaining host metabolic homeostasis, the degree of stability of the core microbiome across a range of environmental conditions emerges as a key trait for predicting the resilience of host populations in aquatic animals [21][22][23] .
In coral reef fishes, recent studies have suggested that intestinal microbiomes influence key physiological functions associated with nutrient acquisition, metabolic homeostasis, and immunity [24][25][26][27] . For example, gut bacteria provide many herbivorous fish hosts with the ability to digest complex algal polymers 24,26,27 and appear susceptible to human disturbances such as eutrophication 28 . The gut microbiome is also a major factor in the innate immune responses to a wide variety of pathogenic microorganisms and other stressors in the surrounding environment 29,30 . Given the rapid physical, chemical, and biotic changes affecting coral reefs, especially in the light of increasing mass coral bleaching events 31 , it is essential to gain a better understanding of how fish gut microbiome assemblages respond to environmental variation so that we can assess how these mutualisms govern host health and resilience to habitat change.
Here, we examined the variability and composition of the gut microbiome of the facultative coral-feeding foureye butterflyfish, Chaetodon capistratus (Linnaeus, 1758), inhabiting a set of reefs that differ markedly in coral cover and diversity across a tropical coastal lagoon (Bahía Almirante) at Bocas del Toro on the Caribbean coast of Panamá. The Chaetodontidae family (butterflyfishes) is among the largest and most iconic families of coral reefassociated fishes 32 and an ideal group for studying microbiome responses to habitat degradation. Chaetodontids range from extreme diet specialists to facultative corallivores and generalists capable of consuming different types of prey such as corals, algae, polychaetes or crustaceans [33][34][35] . Due to their intimate link to the reef benthos, specialized coral-feeding species of Indo-Pacific butterflyfishes are highly sensitive to reductions in coral cover [36][37][38] . Chaetodon capistratus is the only one of the four Western Atlantic Chaetodon species with a relatively high proportion of anthozoans in its diet (mainly hard and soft corals) [39][40][41] . Because of this relative specialization, we chose it as a model species to study relationships between reef habitats and fish host gut microbiomes.
The Bahía Almirante encompasses an inner bay of protected reefs subjected to seasonally high temperatures and a watershed delivering nutrients from agriculture and sewage. In 2010, the bay faced an unprecedented hypoxic event, which led to massive coral bleaching and mortality on some sheltered reefs while others located near the bay's mouth remained unaffected 42 . We capitalized on this gradient of habitat states across the bay to detect variation across fish gut microbiomes in relation to coral degradation. We hypothesized that fish residing on more degraded reefs (i.e., low live coral cover) have a more diverse and variable microbiome as a result of alternative feeding behaviors and potentially increased stress 43 . In contrast, given its role in sustaining host biological functions, we expected that the core microbiome would remain consistent across the habitat gradient.

Results
Benthic habitat and fish density. Reefs located within the three zones classified a priori as outer bay, inner bay, and inner bay disturbed (Fig. 1a), differed in terms of their benthic composition (PCoA; Fig. 2a) with marked differences in the level of live coral cover ( Fig. 2b and Table S1). Live hard coral cover ( Fig. 2b and Table S1) and coral diversity (Shannon diversity; Fig. S1) were highest on reefs of the outer bay. Both stony coral species (i.e., Acropora cervicornis and Agaricia tenuifolia) and fire corals (i.e., Millepora alcicornis, Millepora complanata) dominated at outer bay reefs. In the inner bay zone, reefs displayed an intermediate level of live coral cover ( Fig. 2b and Table S1), largely dominated by the lettuce coral Agaricia tenuifolia. Sponges represented more than a quarter of the benthic cover on these reefs (Fig. S2 and  Table S1). Live coral cover was lowest in the inner bay disturbed zone (Fig. 2b) where dead coral skeleton was prevalent together with sponges ( Fig. S2 and Table S1). Our focal species Chaetodon capistratus showed significantly lower mean density levels at the outer bay than in the two inner bay zones (Fig. S3). Density levels were similar (1-5 individuals per 100 m 2 transect) across all surveyed reefs inside of the bay apart from Cayo Hermanas (SIS, inner bay zone) where up to 25 individuals were recorded in one of the transects (Fig. S3).
Composition of the whole gut microbiome. A total of 5,245,987 high-quality reads were retained for subsequent statistical analyses. The number of reads per sample ranged from 10,369 to 79,466, with a mean ± SD of 41,307 ± 10,990 reads. We identified 10,711 different ASVs in the total dataset. The number of ASVs per sample ranged from 13 to 1,281, with a mean ± SD of 179 ± 210 ASVs. This dataset primarily comprised ASVs belonging to 15 bacterial phyla (Fig. S4a). As predicted, C. capistratus' gut microbiome composition was distinct from the microbiome in seawater and the microbiome of potential prey items (sessile invertebrates) (Fig. S4a, b). Chaetodon capistratus' overall gut microbiome was dominated by Proteobacteria (mainly Gammaproteobacteria, 68.6%) followed by Firmicutes (16.1%), Spirochetes (9.27%), Cyanobacteria (3.98%) (Fig. S4a). Bacteria in the phylum Proteobacteria (Alpha-, Delta-, and Gammaproteobacteria) were dominated by a single genus (Endozoicomonas) in the gut of C. capistratus (93.9%) (Fig. S4b). Endozoicomonas were also abundant in hard-and soft coral samples (23.36 and 41.25% respectively. Firmicutes was abundant in fish guts (16.1% of fish gut bacteria) but representatives of this phylum were nearly absent from potential prey and seawater (Fig. S4a, b). Venn diagrams revealed that fish gut microbiomes shared largely similar proportions of ASVs with coral and sponge microbiomes in each zone (Fig. S5a-c). Fish shared a slightly higher proportion of ASVs with corals (hard and soft coral microbiomes combined) in the inner bay (6.6%; Fig. S5b) and inner bay disturbed zones (6.38%; Fig. S5c) than at the outer bay (5.43%; Fig. S5a). The proportion of shared ASVs between fish gut-and sponge microbiomes was lowest at the inner bay (2.58%; Fig. S5b) and highest at the inner bay disturbed zone (3.41%; Fig. S5c).
Blastn searches against nr/nt NCBI database revealed that ASVs identified as part of the core gut microbiome (i.e., Endozoicomonas) were previously found in scleractinian and soft coral tissue at our study area and in Curaçao among other locations (Table 1). Some Endozoicomonas ASVs were closely related to sequences identified previously in sponges, clams, ascidians, tunicates, and coral mucus as well as the intestinal tract of a coral reef fish species (Pomacanthus sexstriatus). Sequences assigned to Ruminococcaceae closely resembled bacteria reported from herbivorous marine fishes (Kyphosus sydneyanus, Naso tonganus, Acanthurus nigrofuscus, and Siganus canaliculatus), the omnivorous coral reef fish Pomacanthus sexstriatus and a freshwater fish. An Epulopiscium ASV matched to a sequence detected in the guts of two coral reef fishes, the omnivore Naso tonganus and the carnivore Lutjanus bohar and to sequences found in the coral Orbicella faveolata. Other Lachnospiraceae bacteria found in this study resembled sequences known from cattle rumen, hot springs, farm waste, human and other animal feces. Within Ruminococcaceae in Firmicutes, ASVs assigned to the genus Flavonifractor closely resembled bacteria reported from the hindgut of the temperate herbivorous marine fish Kyphosus sydneyanus in New Zealand. Brevinema sequences similar to ours have been previously isolated from the gut of the coral reef fish Naso tonganus as well as freshwater and intertidal fish intestinal tracts. Retrieved Vibrionaceae (genus Vibrio) were similar to sequences found in a coral reef fish gut of Zebrasoma desjardinii. A Romboutsia ASV (family Peptostreptococcaceae), a recently described genus of anaerobic, fermentative bacteria associated with the intestinal tract of animals including humans [45][46][47] , which also occurs in mangrove sediments 48 , matched a sequence found in the tissue of the sea fan Gorgonia ventalina at our study site Bocas del Toro (Table 1).
Alpha diversity of the whole gut microbiome. We estimated alpha diversity using Hill numbers of three different orders of diversity (Hill numbers, {q = 0, 1, 2}) that place more or less weight on the relative abundance of ASVs. This approach allowed for balancing the representation of rare ASVs that might be the result of sequencing errors. Diversity of the whole gut microbiome was lower in fish of the outer bay zone than in fish of the inner bay and inner bay disturbed zones (Fig. 3a-c). Diversity differed significantly among the three zones when taking into account ASV frequency with the Shannon index (Fig. 3b) and when emphasizing abundant ASVs with the Simpson index ( Fig. 3c and Table S3). However, observed ASV richness did not significantly differ among zones ( Fig. 3a and Table S3). Benjamin-Hochberg corrected post hoc tests showed significantly higher Shannon diversity in fish guts of the inner bay zone versus the outer bay zone (Table S4). Fish of the inner bay disturbed zone had a higher microbial diversity than fish of the outer bay zone based on both Shannon and Simpson (Table S4). Pairwise comparisons of alpha diversity between reefs revealed that fish that resided on the reef with the highest level of coral cover overall (37.07%), Salt Creek (SCR, outer bay), had a significantly lower diversity of microbes in their guts than fish from all three inner bay disturbed reefs (RNW, PST, and PBL) for both Shannon and Simpson diversity (Table S5).
Alpha diversity of the core gut microbiome. Diversity of ASVs in the core microbiome was lowest at the outer bay when comparing ASV richness among fish of the outer bay, inner bay, and inner bay disturbed zones and was highest in fish in the inner bay disturbed zone with both the Shannon index and Simpson index ( Fig. 3d-f). Alpha diversity differed significantly among the three zones (Table S3) and pairwise testing revealed that this was largely due to differences between fish of the outer bay and inner bay disturbed zones (Table S4). When compared by reef, lower core microbial diversity in fish from Salt Creek (SCR, outer bay) than fish from other reefs across all zones was responsible for the most significant comparisons (Table S5).
Beta diversity of the whole gut microbiome. Permutational Analysis of Multivariate Dispersion (PERMDISP2) indicated no difference in variability in the whole fish gut microbiome across zones and reefs using dissimilarity metrics that put limited weight on abundant ASVs (Fig. 4a, b and Table S6). However, Bray-Curtis, which more heavily weighs abundant ASVs, identified significantly higher multivariate dispersion for fish from the inner bay disturbed zone than for fish from the outer bay zone ( Fig. 4c and Table S6). The same pattern was observed with phylogenetic dissimilarity metrics. Only the two metrics taking into account relative abundances (i.e., GUniFrac and WUniFrac) revealed significant differences in dispersion patterns among the three zones. Using GUniFrac, an index that adjusts the weight of abundant ASVs based on tree branch lengths, gut microbiomes of fish from the inner bay disturbed zone were significantly more spread in multivariate space than gut microbiomes of fish from the outer bay zone ( Fig. 4e and Table S6). Gut microbial communities were significantly more variable in fish from the inner bay zone than in fish from the outer bay zone using both GUniFrac and WUniFrac (Fig. 4e, f and Table S6).
The three PERMANOVA models explained a small portion of the variance in the composition of the whole gut microbiome using all metrics (2.29-9.22%; Fig. 5a and Table S8). Nevertheless, gut microbiome composition was significantly different between fish from all three zones (zone model), between fish collected inside and outside the bay (position model) and between fish collected on inner bay reefs that differ in coral cover (cover model) when using Jaccard, modified Gower and Bray-Curtis distances ( Fig. 5a and Table S8). Whole gut microbiomes differed using phylogenetic metrics UniFrac and GUniFrac but not when emphasizing microbial relative read abundance (WUniFrac) ( Fig. 5a and Table S8). Pairwise Adonis with Bonferroni corrected P values revealed significant differences among all pairs of zones using non-phylogenetic metrics (Table S10). Pairwise tests were significant using the Unifrac distance except between gut microbiomes of fish from the inner bay and inner bay disturbed zones. None of the pairwise tests using GUnifrac and WUnifrac were significantly different among zones (Table S10).
Whole fish gut microbiomes featured differential relative read abundances across reefs of the inner bay disturbed, inner bay, and outer bay zones (Fig. S7). Gut microbiomes of fish from the inner bay disturbed zone had a lower proportion of microbial reads assigned to Endozoicomonadaceae (Proteobacteria), but a higher proportion of Vibrionaceae and Rhodobacteraceae. In contrast, the relative contribution of Spirochetes and Firmicutes was highest in the guts of fish in the inner bay disturbed zone (Fig. S7). Within Spirochetes, the relative abundance of Brevinemataceae was highest in gut microbiomes of fish from the inner bay disturbed zone, while Clostrideaceae within Firmicutes contributed more to gut microbiomes of fish on inner bay reefs but relatively little to the gut microbiomes of fish in the outer bay zone. Shewanellaceae (phylum Proteobacteria)  represented a higher proportion of the gut microbiome of fish on inner bay disturbed reefs (Fig. S7).
Beta diversity of the core gut microbiome. Patterns in multivariate dispersion were largely consistent between whole and core gut microbiomes. Differences among the three reef zones were significant for metrics that place more weight on ASV relative abundance (common ASVs) (Fig. 4h, i and Table S7). The variability of the core gut microbiome differed significantly between fish from the inner bay and inner bay disturbed zones and between fish from the inner bay disturbed and outer bay zones with highest variability levels in the inner bay disturbed zone. However, none of the phylogenetic metrics showed significant differences in dispersion among zones (Fig. 4j-l and Table S7) As with the whole gut microbiome, the three PERMANOVA models explained a limited amount of variance in the composition of the core gut microbiome (Fig. 5b and Table S9). Yet, composition differed significantly among fish from the three zones (zone model), between fish in the inner bay and outer bay zones (position model) as well as between zones of differential coral cover within the bay (cover model) ( Fig. 5b and Table S9). The core gut microbiome appeared largely similar in composition using all phylogenetic metrics but Unifrac (Table S9). Similar to the whole microbiome, pairwise Adonis tests with Bonferroni corrected P values showed significant differences among almost all pairs of zones when using taxonomic metrics (Table S11). Of the phylogenetic metrics, the only significant differences were found between the inner versus outer bay, and between the inner bay disturbed versus outer bay zone using Unifrac (Table S11). Differences in the composition of the core microbiome among reef zones was largely driven by changes in the relative abundance of ASVs assigned to the genus Endozoicomonas (class Gammaproteobacteria) (Fig. S5b). For example, the most common Endozoicomonas ASV (ASV1) was much more represented in the guts of fish from the outer bay and inner bay zones than in the guts of fish in the inner bay disturbed zone, while Endozoicomonas ASVs relative abundances appeared more evenly distributed towards the inner bay disturbed zone. In contrast, bacteria in the genus Brevinema (phylum Spirochetes) were most abundant relative to other members of the core in fish of the inner bay disturbed zone and least abundant in the outer bay zones. The giant bacterium Epulopiscium (family Lachnospiraceae, order Clostridia), which is known to aid the digestion of algae in surgeonfishes, contributed more to the core gut microbiome of fish on reefs in the inner bay disturbed zone than the inner and outer bay zones (Fig. S5b).
Prevalent ASVs in each reef zone. A machine learning-based, de-noising algorithm (PIME) was used to detect sets of ASVs in the whole gut microbiome that significantly contribute to differences between reef zones. The initial out-of-bag (OOB) error rate (i.e., the prediction error in a Random Forest model) for our unfiltered dataset was greater than 0.1 (PIME, OOB 0.27) indicating that PIME filtering would effectively remove noise. PIME identified a prevalence cut-off of 65% for the highest improved accuracy (OOB = 2.25) indicating that the model was 97.75% accurate (Table S12). The validation step showed that randomized errors (Fig. S8b) corresponded with the predicted prevalence cut-off value of 0.65 indicating the absence of false positives (Type I error).
After selecting ASVs that were present in at least 65% of the fish guts within each zone, the filtered dataset comprised 17 ASVs in eight families; i.e., Endozoicomonadaceae, Ruminococcaceae, Pirellulaceae, Lachnospiraceae, Brevinemataceae, Cyanobiaceae, Rhodobacteraceae, and Peptostreptococcaceae ( Fig. 6 and Tables S12, S13, S14). Fish of the inner bay zone showed the highest richness levels with 13 ASVs, compared to eight and nine ASVs in fish of outer bay and inner bay disturbed zones, respectively (Fig. 6). An Endozoicomonas ASV (ASV1), which was also a dominant component of the core, had a much higher relative abundance in fish of the outer bay zone than in fish of the inner bay disturbed zone (Fig. 6). Communities differed most in composition between fish of the outer bay and inner bay disturbed zone, whereas fish of the inner bay zone reflected an intermediate community between these two comprising the highest richness of Endozoicomonas with six ASVs. As in the core community, the Endozoicomonas assemblage was slightly less diverse in the disturbed zone (three ASVs) and featured more similar relative read abundances than in the outer bay zone (four ASVs) where a single ASV was dominant. Two distinct ASVs of the giant bacterium Epulopiscium (family Lachnospiraceae) were prevalent in fish in both the inner and inner bay disturbed zones but were more abundant on disturbed reefs. Disturbed reefs uniquely featured anaerobic gut bacteria in the genus Romboutsia (family Peptostreptococcaceae) (Fig. 6).

Discussion
Detecting how the spatial turnover of microbiomes varies within and among host populations, and in relation to habitat characteristics is essential to understanding and predicting the response of host species to environmental change. We show that whole gut microbial communities were significantly more diverse and variable in fish from inner bay disturbed reefs than from the outer bay zone in terms of some but not all measured diversity components. Conspicuously, the core microbiome, a small set of microbial strains that may form sustained relationships with the fish host, also showed higher dispersion on degraded reefs suggesting greater variability of microbial assemblages among individual fish. Significant differences in diversity and group dispersion were observed mostly in the relative abundances of the frequent and common taxa. Highly variable host-associated microbial communities have been observed in humans with immunodeficiency syndromes (reviewed in ref. 49 ) and in marine animals such as scleractinian corals and anemones under acute stress 43,50-52 . Zaneveld et al. 43 referred to this pattern of variability as the "Anna Karenina principle" applied to host-associated microbiomes (AKP). They argued that this is a common but often overlooked response of organisms that become unable to regulate  Reductions in coral cover may increase foraging costs if, for example, fish spent more energy to search, capture, and handle their prey. Indeed, physiological stress imposed by environmental conditions may cause immune signals that imbalance the gut microbiome 30,43,53,54 . Disturbance to the microbiome, in turn, can affect the brain and further alter behaviors related to movement such as the ability to forage 30,55 . The low coral cover may also increase stress through intra-and interspecific competition for resources. For example, social stress in the form of aggressive interactions among conspecifics was shown to alter the behavior and microbial assemblages associated with mice by setting off immune responses critical to host health 56 . In Indo-Pacific obligate feeding butterflyfishes coral degradation was shown to decrease aggressive encounters among and within Chaetodon species 57 as well as change the frequency of pair formation 58 , and the way species responded to the loss of the coral resource depended on their level of dietary specialization 57,58 . Foraging on degraded reefs may also increase predation risk when architectural complexity is reduced 59 . Anxiety-like behaviors induced by exposure to predators can lead to sustained physiological stress in vertebrates (reviewed in ref. 60 ).
Another possible explanation for more variable gut microbiomes on disturbed reefs could be increased behavioral heterogeneity among fish individuals (e.g., feeding behavior). Where preferred food sources are scarce, foraging behavior may become more diverse and lead to increased individual specialization on alternative food items 61,62 translating into more varied gut microbiomes. In this scenario, the higher variation in gut a Whole microbiome b Core microbiome  Fig. 6 PIME filtering zones 65% prevalence. Comparison of fish gut microbiomes among three reef zones. The whole fish gut microbial dataset was filtered using Prevalence Interval for Microbiome Evaluation (PIME) 134 to detect which ASVs were responsible for differences among zones. Using machine learning, PIME de-noises the data by reducing within-group variability. Based on the algorithm, we selected a 65% prevalence cut-off resulting in a filtered dataset of 17 ASVs at a low error rate (OOB = 2.25) and high model accuracy (97.75%). microbial assemblages would be the result of behavioral adjustments (acclimatization) to alternative habitat conditions without necessarily causing stress. Higher alpha diversity across fish gut microbiomes in the inner bay disturbed zone supports this explanation. Although C. capistratus is able to consume a broad range of diet items 39,63,64 , deviations from its preferred coral prey may come with fitness consequences as shown for Indo-Pacific butterflyfishes 33,65,66 . For example, other authors found that obligate corallivorous Chaetodon species have reduced energy reserves at reefs where they diversify or shift their diet in response to limited coral availability 33,66 .
Apart from patterns of microbiome variability, the significant differences in the composition of the whole gut microbiome (as opposed to the core microbiome) in nearly all comparisons (i.e., among all three zones, between inner and outer bay, and between inner bay disturbed and undisturbed) may primarily reflect changes in diet. Specifically, in the inner bay disturbed zone where coral cover was low, the microbial assemblage suggests (i) potential changes in the invertebrate prey community and (ii) a more broad, likely omnivorous trophic profile indicated by a distinct Endozoicomonas community in codominance with anaerobic fermentative bacteria. The increased prevalence of fermentative microbes at disturbed reefs might reflect the consumption of algae and potentially sponges given their high availability in this zone. However, we lack information on the effects of sponge consumption on fish gut microbiomes. Epulopiscium, often considered a host-specific symbiont of herbivorous surgeonfishes (family: Acanthuridae) 26,27,67,68 , was present in the core microbiome and identified as distinct to the inner bay with predominance at disturbed reefs. This may suggest that C. capistratus can assimilate nutrients from algae and that this metabolic function is enhanced on degraded reefs by the increase in key microbial functional groups. Alternatively, the fish in our study may take up these microbes while foraging for invertebrates on the epilithic algal matrix. Overall, levels of Epulopiscium were approximately similar to those previously found in omnivores and detritivores in the Red Sea 26 with the two most prevalent ASVs matching a strain previously isolated from the turf algal grazer Naso tonganus 69 . Additionally, the presence of Rhodobacteraceae, which are often found associated with algal biofilms 70,71 , may suggest detritivory but might also be related to the consumption of mucus from stressed 72 and diseased corals where Rhodobacteraceae are also found 73,74 . The lower relative abundance of a compositionally distinct Endozoicomonas community on disturbed reefs could reflect different proportions of prey species featuring Endozoicomonas 75,76 in the diet of C. capistratus.
In contrast, a single dominant Endozoicomonas ASV along with a few Firmicutes characterized the gut microbiome of C. capistratus on outer bay reefs. The presence of some Endozoicomonas ASVs shared between fish guts and potential prey (i.e., hard corals, soft corals, zoanthids, sponges), including matches to microbial sequences previously detected in two coral species (Orbicella faveolata and Poritis asteroides) in our study area at Bocas del Toro 73,77 , suggests the horizontal acquisition of these microbes via feeding on corals. In addition, we identified an ASV in the genus Ruegeria as indicative of the outer bay reefs, which matched a sequence previously retrieved from the soft coral species Pterogorgia anceps on the Caribbean coast of Panamá (GenBank Accession: MG099582) and which was also present across samples of potential prey taxa including hard and soft corals and sponge-infauna. Endozoicomonas originating from the food could potentially play a role in promoting the assimilation of nutrients via interactions with resident bacteria 78 .
The core microbiome composition significantly differed across the inner bay between fish from disturbed and undisturbed reefs where environmental conditions are generally homogeneous except for the proportion of live coral cover. Despite our models accounting for relatively little variance, this finding may suggest that bacterial communities that are most likely to have intimate metabolic interactions with C. capistratus might fail to provide important functions to hosts in severely degraded habitats. However, we cannot exclude the contribution of other factors that were not measured in this study such as microbial plasticity mediated by diet, gut colonization history 79 , and/or potential genetic differentiation between the inner bay and outer bay fish populations [80][81][82] .
Our analysis identified ten Endozoicomonas ASVs as part of the core microbiome indicating potential true resident symbionts. Members of the genus Endozoicomonas spp. are known as bacterial symbionts of marine sessile and some mobile invertebrates and fishes 44,76,83-87 . Reverter et al. (2017) 86 found Endozoicomonas associated with butterflyfish gill mucus in Chaetodon lunulatus and Parris et al. (2016) 87 found Endozoicomonas in the gut of damselfishes (family: Pomacentridae) and cardinalfishes (family: Apogonidae) pre-and to a lesser extent post-settlement on the reef. Corallivory in butterflyfishes has evolved in close association with coral reefs 32,88 and this likely involved adaptive mechanisms to metabolize defense compounds from corals and many other sessile invertebrates (e.g., polychaetes). Adapted gut microbial communities may help butterflyfish hosts cope with toxins or facilitate the digestion of complex prey tissues as in insects 89 , mammalian herbivores 90 , and surgeonfishes 26 . It is likely that the gut microbial profile of C. capistratus -featuring high abundance Endozoicomonas-facilitates the digestion of complex coral prey. More detailed knowledge will be required to understand whether the potential intake of Endozoicomonas via fish browsing on sessile invertebrates plays a role in nutrient uptake in trophic strategies such as fish corallivory.
We detected increases in gut microbiome variability, diversity, and spatial community turnover. These patterns extended to the core microbiome suggesting signs of potentially altered functioning that may affect fish hosts on reefs with extremely low levels of live coral cover. Nonetheless, the density of C. capistratus was comparable across both inner bay zones indicating that the lack of live coral cover may not immediately impact the persistence of populations. Significantly lower densities at the outer bay may potentially relate to spatial patterns of larval recruitment 91 and/or differences in wave exposure across reefs affecting the energy expenditure fish allocate towards swimming performance and feeding 92 . Additional work should focus on linking changes in the gut microbiome to direct measures of diet and host health. Our results give insight into the poorly understood spatial fluctuations in host-associated microbial communities across a natural system. This work highlights intricate links between ecosystem-scale and microbial-scale processes, which have so far been mostly overlooked. We suggest there is an urgent need to integrate measurements of the role of microbes in the response of reef fishes to the global loss of coral reefs.

Methods
Study area. Bahía Almirante, located in the Bocas del Toro Archipelago on the Caribbean coast of Panamá, is a coastal lagoonal system of~450 km 2 where numerous, relatively small, and patchy fringing coral reefs occur 93 . Hydrographic and environmental conditions vary across the semi-enclosed bay but are generally characterized by limited water exchange with the open ocean 42 . Furthermore, areas of the bay are subjected to uncontrolled sewage and dredging due to increasing coastal development and agricultural runoffs from the adjacent mainland [94][95][96][97] . A total of nine discontiguous reefs distributed from the mouth towards the inner bay were selected for this study based on distinct hydro-geographical zones and disturbance history, resulting in three distinct reef zones with three replicates each (total n = 9 reefs) (Fig. 1a). Throughout the manuscript, we will refer to these three distinct reef zones as "outer bay", "inner bay", and "inner bay disturbed". Outer bay reefs [Salt Creek (SCR), Cayo Corales (CCR), and Popa (PPR)] are located at the mouth of the bay marking a transition zone between the inner bay and the open ocean. These reefs represent typical Caribbean reef communities featuring both massive and branching coral colonies with higher benthic cover and diversity as compared to the inner bay (Fig. 1b). Inner bay reefs [Almirante (ALR), Cayo Hermanas (SIS), and Cayo Roldan (ROL)] are largely coral and sponge dominated reefs and have lower coral diversity than the outer bay reefs (Fig. 1c). Inner bay disturbed reefs [Punta Puebla (PBL), Punta STRI (PST) and Runway (RNW)] were heavily impacted by the 2010 hypoxic event 42 , which resulted in the current cover of largely dead coral comprised of formally prevalent Agaricia and Porites species (Fig. 1d). Prior to this disturbance, both study zones located inside of the bay exhibited comparable benthic communities of similar health states. For example, the Punta STRI reef (PST) at the now disturbed zone featured 26.9% coral cover in 2005 98 .
Benthic habitat and fish density. Visual surveys of benthic cover and focal fish species densitiy were conducted between May and June 2016. At each of the nine reefs, three 20 m transects were placed parallel to the shore at 2-4 m depth. Benthic community cover was estimated from 100 cm × 70 cm photographic quadrats placed every 2 m, resulting in a total of 10 quadrats per transect. Photos were analyzed on the CoralNet platform 99 using a stratified random sampling design (10 rows × 10 columns with 1 point per cell for a total of 100 points per image). The first 15 of all photos were manually scored to train the algorithm. The remaining photos were then processed by the automated assignment tool, and assignments were subsequently verified for each point. Due to the difficulty involved with photo-based taxonomic identification, analyses were conducted at the level of broad benthic categories which comprised the following: live hard coral, dead hard coral, live soft coral, sponge, zoanthid, other invertebrates, seagrass, grazable substrate, macroalgae, rubble, sand, shade and "unknown". Within the live hard coral, dead hard coral, and live soft coral categories, identification was done at the genus or species level when possible. The mean cover of each benthic category was calculated per reef. To estimate focal fish species density, C. capistratus individuals were counted along each 20 × 5 m belt transects used subsequently for the benthic surveys while swimming slowly using scuba (except at CCR).
Sample collection. The foureye butterflyfish, Chaetodon capistratus, is a common member of Caribbean coral reef fish communities (IUCN classified as least concern) 100 with a distribution that extends across the subtropical Western Atlantic 101,102 (Fig. 1e). The following protocol of fish capture and euthanization had been approved by the Smithsonian Tropical Research Institute's Institutional Animal Care and Use Committee (IACUC). An average of 11 individual adult fish were collected at each of the nine reefs (min = 7; max = 16; total = 102) by spearfishing in February and March 2018 (Table S15). Captured fish were immediately brought to the boat, anesthetized with clove oil, and placed on ice in an individual and labeled sterile Whirl-Pak bag. Upon return to the research station, fish were weighed (g wet weight), and both standard length (mm SL) and total length (mm TL) were measured with a digital caliper. The intestinal tract of each fish was removed under a laminar flow hood using tools decontaminated with 10% sodium hypochlorite. The intestinal tracts were then preserved in 96% ethanol in individual 15 ml or 5 ml centrifuge tubes and stored at −20°C until DNA extraction. To assess microbial communities present in the fish's environment, we also obtained samples of seawater and potential prey taxa. At each of nine reefs, a total of four liters (2 × 2 L at each reef) of seawater was collected immediately above the reef substratum using sterile Whirl-Pak bags and filtered through a 0.22 µm nitrocellulose membrane (Millipore) and a total of 18 seawater samples was included into downstream analysis. Small pieces of hard coral (Siderastrea siderea, n = 2; Porites furcata, n = 2; Agaricia tenuifolia n = 2), soft coral (Antillogorgia bipinnata, n = 1; Plexauridae sp.; n = 1), sponges (Amphimedon compressa, n = 1; Chondrilla caribensis, n = 1; Demospongiae spp., n = 4), macroalgae (Amphiroa sp., n = 1), turf (n = 1), and zoantharia (Zoanthus pulchellus, n = 1; Palythoa caribaoerum, n = 1) were collected opportunistically at least at one of the three habitat zones and kept in sterile Whirl-Pak bags on ice on the boat. At the field station, samples were individually placed in 50 ml or 15 ml centrifuge tubes with 96% ethanol and stored at −20°C until DNA extraction.
DNA analysis. The mid-and hindgut of the gastrointestinal tract of each fish was opened longitudinally to isolate the digesta and the mucosa tissue by lightly scraping the intestinal epithelium. Between 0.05 and 0.25 g of both tissue types combined was used for DNA extraction using the Qiagen Powersoil DNA isolation kit following the manufacturer's instructions with minor modifications to increase yield (see supplementary methods). Each tissue sample of potential prey organisms (invertebrates and macroalgae) was homogenized in separate vials. Additionally, infaunal communities (small worms) were isolated from two sponges, Amphimedon compressa and Dysidea sp. and the tissue homogenized for each sponge separately. DNA was extracted (0.25 g per sample) following the same protocol as described for the intestinal microbiomes. Seawater DNA was isolated from nitrocellulose membrane filters using the Qiagen Powersoil Kit following a modified protocol described in ref. 103 .
A dual Polymerase Chain Reaction (PCR) approach was used to amplify the V4 hypervariable region (primers 515F 104 and 806R 105 ) of the 16 S ribosomal rRNA gene of each sample (Table S16). Subsequently, the product of all samples was sequenced by combining them into a single Illumina MiSeq sequencing run. Our protocol followed the 16 S Illumina Amplicon Protocol of the Earth Microbiome Project 106 using locus-specific primers to which Illumina "overhang" sequences were appended. These overhang sequences served as a template to add dual index Illumina sequencing adapters in a second PCR reaction (see supplementary methods for detailed PCR protocols). The final product was sequenced on the Illumina MiSeq sequencer (reagent kit version 2, 500 cycles) at the Smithsonian Tropical Research Institute with a 20% PhiX spike. The absence of contaminants was confirmed with negative DNA extractions and negative PCR amplifications (see supplementary methods for detailed DNA extraction and PCR protocols).
Analysis of sequence data. All analyses were conducted with the statistical software R version 3.6.158 107 . Illumina adapter and primer sequences were removed from forward and reverse reads using Cutadapt 108 with a maximum error rate of 0.12 (-e 0.12). The remaining reads were filtered and trimmed based on their quality profiles and potential chimeras were removed using DADA2 version 1.12.1 109 . Sequences were discarded if they had more than two expected errors (maxEE = 2), at least one ambiguous nucleotide (maxN = 0), or at least one base with a high probability of erroneous assignment (truncQ = 2). Forward and reverse reads were trimmed to 220 and 180 bp respectively to remove lower quality bases while maintaining sufficient overlap between paired-end reads. Sequences were kept when both the forward and reverse reads of a pair passed the filter. Quality filtered reads were de-replicated and Amplicon Sequence Variants (ASVs) inferred. Paired-end reads were merged and pairs of reads that did not match exactly were discarded. Taxonomy was assigned to each ASV using a DADA2 implementation of the naive Bayesian RDP classifier 110 against the Silva reference database version 132 111 . ASVs identified as chloroplast, mitochondria, Eukaryota, or those that remained unidentified at the kingdom level were removed from the dataset. Sequences of each ASV were aligned using the DECIPHER package version 2.0 112 . The PHANGORN package version 2.5.5 113 was then used to construct a maximum likelihood phylogenetic tree (GTR + G + I model) using a neighbor-joining tree as a starting point. Fourteen samples containing few sequences (<10,000) were removed from the dataset (Fig. S9). The remaining samples were rarified without replacement to even sequencing depth (n = 10,369 sequences) for downstream analysis. Our approach followed the recommendation for the normalization of sequencing data 114 . Statistical analysis was conducted using phyloseq version 1.28.0 115 .
Delineation of the core gut microbiome. To identify the persistent bacteria associated with the fish gut (i.e., the "core microbiome" 19,116 ) including taxa that might be potentially beneficial to the fish host, we employed a statistical approach taking into account both relative abundance and relative frequency of occurrence of ASVs as opposed to the common procedure of using an arbitrary minimum frequency threshold based on presence-absence data only 116 . Indicator species analysis 117 (labdsv package version 2.0-1) 118 was used to identify which ASVs were relatively more abundant and predominantly found in fish guts and not in their surrounding environment. We calculated an Indicator Value (IndVal) Index between each ASV and two groups of samples: (1) all fish gut samples, and (2) all seawater and sessile invertebrate samples upon which fish potentially feed. Statistical significance of the association between ASVs and groups of samples was tested using 1000 permutations. ASVs were considered indicators of fish guts (i.e., components of the core) based on a maximum probability of P value = 0.01. Sequences of ASVs identified as part of the core microbiome were compared to the non-redundant nucleotide (nr/nt) collection database of the National Centre for Biotechnology Information (NCBI) using the Basic Local Alignment Search Tool for nucleotides (BLASTn) 119 . We extracted metadata associated with all sequences that matched each query at 100% similarity or the first five top hits to identify within what taxa, environment, and/or habitat each core ASV and close relatives had been previously found.
Diversity analysis. The workflow of our microbial community analysis is visualized in a diagram (Fig. 7). To account for the presence of rare sequence variants caused by sequencing errors or other technical artifacts 120 , we used Hill numbers 121 following the approach recommended by ref. 122 for sequence data to compare alpha diversity among groups of samples. Hill numbers allow scaling the weight put on rare versus abundant sequence variants while providing intuitive comparisons of diversity levels using "effective number of ASVs" as a measuring unit [121][122][123] . This approach allowed for balancing the over-representation of rare ASVs that might be inflated due to sequencing errors 124 . We calculated three metrics that put more or less weight on common species: (1) observed richness, (2) Shannon exponential that weighs ASVs by their frequency, and (3) Simpson multiplicative inverse that overweighs abundant ASVs. Alpha diversity was calculated and visualized using boxplots for the whole and core fish gut microbiomes. Because Shapiro-Wilk tests indicated that the data were not normally distributed, non-parametric Kruskal-Wallis tests were used to compare alpha diversity among reefs (n = 9) and the three reef zones (outer bay, inner bay, and inner bay disturbed) with post hoc Dunn tests.
To test the hypothesis that fish gut microbiomes are more variable between individuals at disturbed sites, we calculated non-parametric Permutational Analysis of Multivariate Dispersion (PERMDISP2) (betadisper function, vegan package implemented in phyloseq 115,125 ). PERMDISP2 is a measure of the homogeneity of variance among groups and compares the average within-group distance to centroid between each predefined group of samples in multidimensional space. We used a range of phylogenetic and non-phylogenetic dissimilarity metrics that differentially weigh the relative abundance of ASVs to identify the effect of abundant ASVs [Phylogenetic: Unifrac, Generalized Unifrac (package GUniFrac) 126 and Weighted Unifrac 127 ; non-phylogenetic: Jaccard 128 , modified Gower with log base 10 129 and Bray-Curtis 130 ]. P values were obtained by permuting model residuals of an ANOVA (Analysis of Variance) null-model 1000 times (betadisper function, vegan 115,125 ).
Differences in microbial composition were tested using Permutational Multivariate Analysis of Variance (PERMANOVA) with the Adonis function in vegan 131 computed with 10,000 permutations. Comparisons were made (1) among fish gut microbiomes of the three reef zones ("zone model"), (2) between fish gut microbiomes of outer bay reefs versus inner bay reefs ("position model") and (3) between fish gut microbiomes of inner bay reefs and inner bay disturbed reefs which differed in coral cover ("cover model"). PERMANOVA is robust to the effect of heterogeneity of multivariate dispersions in balanced or near balanced designs as in our study 132 . Pairwise Adonis with Bonferroni corrected P values was computed using the pairwise Adonis function (version 0.4) 133 .
Finally, we used the Prevalence Interval for Microbiome Evaluation (PIME) package (version 0.1.0) 134 to identify sets of ASVs that are frequently found in fish guts in each zone at the Bahía Almirante (outer bay, inner bay, inner bay disturbed). This analysis is aimed at identifying microbial ASVs that are differentially prevalent among zones. PIME uses a supervised machine learning Random Forest algorithm to reduce within-group variability by excluding lowfrequency sequences potentially confounding community comparisons of microbiome data 134 . PIME identifies the best model to predict community differences between groups by defining a prevalence threshold that retains as many ASVs as possible in the resulting filtered communities (i.e., the random forest classifications) while minimizing prediction error (out-of-bag error, OOB). To do so, the algorithm uses bootstrap aggregating (100 iterations) of each sample group at each filtering step (prevalence interval) by 5% increments. Random Forest calculates a global prediction from a multitude of decision trees based on the bootstrap aggregations and estimates the out-of-bag error rate (OOB) from omitted subsamples during aggregating 135 . Validation was done by randomizing the original dataset (100 permutations) and subsequently estimating Random Forest error to determine if group differences in the filtered dataset were due to chance (pime.error.prediction function, PIME) 134 . A second function (pime.oob.replicate, PIME) 134 repeated the Random Forest analysis using the filtered dataset for each prevalence interval without randomizing group identity. In a preliminary step, we assessed whether the OOB error for our unfiltered data was >0.1, which indicated that de-noising with PIME would improve model accuracy.
Statistics and reproducibility. Statistical analyses were carried out for the 16 S sequencing data and in-situ transect data (benthic photographic quadrats and visual censuses of fish communities). Details allowing the reproducibility of all analyses are provided in the methods section (including sampling sizes and numbers of replicates). A diagram illustrating the statistical workflow is also included (Fig. 7) and the R code is available on the project website: https://github.com/bocasbiome/web/.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Sequencing data has been submitted to the NCBI Short Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra) under bioproject number Accession: PRJNA718434 ID: 718434. Raw data are available on Dryad Digital Repository https://doi.org/10.5061/ dryad.m905qfv28 153 . Fig. 7 Microbial community analysis. Microbial community analysis workflow illustrating how we subsetted the whole fish gut microbiome dataset to delineate the core gut microbiome and gut microbial communities by zone, respectively. To identify the core microbiome, we used indicator analysis 117 between the whole fish gut microbiome and the environmental sample fraction consisting of samples of potential fish prey taxa and the surrounding seawater. Diversity analysis was done for the whole and core fish gut microbiome, respectively. The whole fish gut microbiome was filtered for prevalence with a machine learning-based algorithm (PIME) 134 to detect community differences among zones that reflect fish-microbiome responses to the habitat gradient. Created with BioRender.com. The fish icon is adapted from a color photograph of Chaetodon capistratus obtained from https://biogeodb.stri.si. edu/caribbean/en/pages with permission by D R Robertson. Icons of benthic organisms obtained from the IAN Symbol Libraries: Tracey Saxby and Joanna Woerner, Integration and Application Network (ian.umces.edu/media-library). https://creativecommons.org/licenses/by-sa/4.0/.